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Abstract 

We present density functional theory calculations of phosphorus dopants in bulk silicon and 
of several properties relating to their use as spin qubits for quantum computation. In contrast 
to previous ab initio calculations, we have used an explicit treatment for the phosphorus donor 
and examined the detailed electronic structure of the system as a function of the isotropic doping 
fraction, including lattice relaxation due to the presence of the impurity. Doping electron densities 
{Pdoped — Pbulk) and Spin densities {p^ — p^) are examined in order to study the properties of the 
dopant electron as a function of the isotropic doping fraction. Doping potentials {Vdoped — Vbuik) are 
also calculated for use in calculations of the scattering cross-sections of the phosphorus dopants, 
which are important in the understanding of electrically detected magnetic resonance experiments. 
We find that the electron density around the dopant leads to non-spherical features in the doping 
potentials, such as trigonal lobes in the (001) plane at energy scales of +12 eV near the nucleus and 
of -700 meV extending away from the dopants. These features are generally neglected in effective 
mass theory and will affect the coupling between the donor electron and the phosphorus nucleus. 
Our density functional calculations reveal detail in the densities and potentials of the dopants 
which are not evident in calculations that do not include explicit treatment of the phosphorus 
donor atom and relaxation of the crystal lattice. These details can also be used to parameterize 
tight-binding models for simulation of large-scale devices. 



*This author's current location is AX Division, Lawrence Livermore National Laboratory 
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I. INTRODUCTION 



Dopants in silicon show potential as qubits for solid-state quantum computers l|-|3|], with 
the advantages of scalability as well as the promise of utilizing the existing semiconductor 



industry and its processing techniques jl-5|. 



The theory of Group V dopants such as phos- 



phorus in silicon [6| is useful for describing the quantum nature of the electrons in these 
systems, as well as for developing schemes to circumvent one of the most__cliallenging aspects 

14| . In order to 



of solid-state quantum computers, namely environmental decoherence 
provide a benchmark for such theories and also to use as a starting point for building efficient 
and accurate tight binding methods, an ab initio description of dopants in silicon is desired. 
However, the size of the systems required to describe doped silicon at or near the single- 
dopant limit is large, making such a description computationally expensive. In this work, 
we present large-scale density functional theory calculations for phosphorus-doped silicon 



15 



supercells with up to 432 atoms. We make comparisons to other theoretical works 
to determine what can and cannot be captured by approximate or single-electron theories 
for the doped-silicon systems. 

Previous efforts to describe the electronic structure of silicon dopants include effective 
mass approaches beginning with the work of Kohn and Luttinger |6|] and continuing with 
many others 17|, Il9l-l24j. including Fan g et . al. who perform two-electron Hartree-Fock 
calculations within effective mass theor y |25l| . These efforts include calculating the effects of 
applied electric and magnetic fields 20|, |22[ and the coupling of two donors via the exchange 



interaction 



17 



21 



26 



m- 



, |23| . Tight-binding calculations have also been performed 
eluding a calculation of the quadratic Stark coefficient of the hyperfine interaction which has 
reproduced experimentally measured values more accurately than effective mass theory 



29|. 



Two-dimensional layers of dopants in silicon known as 6 



avers have been described using 



density functional theory with mixed pseudopotentials [16|, |30|, |3l| , which treat the dopant 
and silicon atoms in the layer using the same core potential. These DFT calculations and 
a number of additional calculations (see Ref. 32[ and Refs. [14-24] therein) show a large 
amount of disagreement for calculated properties such as the valley splitting. Additionally, 
electrically detected magnetic resonance (EDMR) 33| has called into question a theoretical 
picture of the scattering of electrons in a two-dimensional electron gas (2DEG) from dopants 



in silicon 



34 



35|. An ab initio description of a dopant in silicon is therefore useful both as a 
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benchmark and for determining the details of the electronic structure of an isolated dopant 
which can subsequently be used to calculate more accurate spin-dependent scattering cross 
sections. Although these are expensive calculations, we have been able to perform large-scale 
calculations using the computational resources at Lawrence Livermore National Laboratory. 



II. METHODS 



We have used the Quantum ESPRESSO suite of programs [36| in order to perform 
density functional theory calculations using a basis set of plane waves. Face-centered cubic 
lattices of substitutionally doped silicon were prepared at variable dopant ratios by substi- 
tuting 1 phosphorus atom in unit cells with 53, 127, 249, and 431 silicon atoms, respectively. 



We used the Perdew-Burke-Ernzerhof (PBE) density functional 



pseudopotential for phosphorus (P.pbe-n-van.UPF from Ref. 39|) and a norm-conserving 



with an ultrasoft 



pseudopotential that was calculated using FHI98PP 40[ for silicon. A plane-wave energy 



cutoff of 65 Ry. was chosen based on convergence of the total energy and pressure in the 



smallest supercell. K-space sampling was performed using a Monkhorst-Pack [4l| grid of 
8x8x8 (54 atom), 6x6x6 (128 atom), 4x4x4 (250 atom), and 2x2x2 (432 atom) grid 
points. As a reference for the energy and properties of bulk silicon, a two-atom silicon cell 
was used with a grid of 20x20x20 k-points. In this cell, one silicon is placed at the origin 
and another is placed at the point (|, |, |), where a is the lattice constant. The basis vec- 
tors of the FCC cell are (|, |, 0), (|, 0, |), and (0, |, |); the two sihcon atoms are repeated 
at every integer multiple of the basis vectors. The lattice constant of the doped supercells 
was determined as multiples of the 5.46 A lattice constant for bulk silicon computed with 
the pseudopotential used in this study. Geometric minimization of the total energy was 
performed for all doped systems studied here. In the calculations with < 432 atoms, all 
of the atoms were allowed to move during the simulation, while for A^ = 432, the atoms 
at the edges of the supercell were frozen in order to make better comparisons to the bulk 
silicon system. The volume of each supercell was constant during the geometric optimiza- 
tion. These computations required 600 processors on the Ansel supercomputer at Lawrence 
Livermore National Laboratory, a 324 node system using the Intel Xeon architecture rated 



at a peak of 43.5 TFLOP/s 42|. A single point energy calculation for the 432-atom cell 



took about 19 CPU hours to complete. 



3 



In Section imZl we use the doping potentials as a measure of the electronic environment 
of the dopants. The electric potential for each cell is obtained by adding the nuclear term 
{Vnuc), coulombic Hartree term {VHartree), and PBE exchange-correlation term (V^c) at the 
converged electronic density Pconv, 

^ ^nuc ~l~ ^Hartree[Pccmv\ ~l~ ^^c[Pcorn;]- (l) 

The density Pconv is that which minimizes the energy of the system according to Eq. ([T]). 
The doping potential is then obtained by subtracting the potential obtained for the undoped 
cell from that of the doped cell, 

^doping ^doped ^undoped' (2) 

Similarly, the doping density discussed in Section [Tll Cl is obtained by subtracting Pconv of the 
undoped cell from that of the doped cell. Finally, the spin densities discussed in Section [III Bl 
are obtained by subtracting the spin-up and spin-down portions of the density Pconv of a 
single cell. 

In order to estimate the exchange coupling (J), we use the DFT broken spin symmetry 



states 



43 



47| 



^ = Ebs - Ehs, (3) 

where BS and HS denote, respectively, the broken symmetry ground state and high-spin 
state in which the spin density is constrained to give a spin of 5* = | in each unit cell. This 
procedure empirically has given exchange couplings with a satisfactory degree of accuracy 



for molecular systems 



43 



44| , and the exchange couplings it calculates can be thought of as 



the spin-coupling parameter J in both Ising and Heisenberg models of the spin interactions 



in the periodic system 48l. |49|. 



III. RESULTS AND DISCUSSION 



Is 



The Bohr radius of phosphorus in silicon is estimated to be about 2.5 nm (6|. In order 
to reach the single dopant limit, the phosphorus atoms in a silicon matrix would have to be 
spaced a minimum of ~ 5 nm apart. This would require a supercell of tens of thousands 
of atoms and is outside the realm of feasibility for most DFT calculations. Here, we have 
therefore assessed the energies and properties of doped silicon as a function of the doping 



fraction approaching this hmit, in order to develop a better understanding of how exphcitly 
treating the phosphorus donor atom and lattice relaxation effects the electronic structure. 
In Table [U we show the convergence of the energy per atom as a function of the size of the 
supercell. We have also compiled the doping energies of each cell relative to that of the 

TABLE I: Convergence of the energy as a function of unit cell size. For each unit cell from 54 to 
432 atoms, we give the energy per atom (^^^, column 2) and the doping energies (columns 3, 5) 
and geometry optimization energies (column 4). Subscripts denote whether a cell is doped (with 1 
phosphorus atom) or undoped (all silicon atoms), and superscripts denote whether the geometry 
is that of bulk silicon or the optimized geometry of the doped cell. 



Unit cell 


^doped 
Natom 

(Rydberg) 


rpbulk rpbulk 

doped undoped 

(Rydberg) 


popt _ pbulk 
doped doped 

(mRydberg) 


popt _ rpbulk 
doped undoped 

(Rydberg) 


Si2 


-7.8814 








SissP 


-8.0668 


-10.0116 


-2.2 


-10.0138 


Sil27P 


-7.9596 


-10.0137 


-2.0 


-10.0158 


81249? 


-7.9214 


-10.0141 


-2.4 


-10.0165 


Si43lP 


-7.9045 


-10.0175 


-2.4 


-10.0198 



equivalent bulk silicon cell, with the doped cell energy calculated for both the bulk silicon 
geometry (column 3) and the optimized geometry of the doped cell (column 5). We will 
refer to the difference in energy between the geometrically relaxed doped cell and the bulk 
geometry undoped cell as the "relaxing and doping energy." This quantity is listed in column 
5 of Table [B The geometric relaxation energy of the doped system (column 4 of Table [T]) is 
similar (around 2 mRy) for each of the cells. The relaxing and doping energy increases in 
magnitude by 6.0 mRy (81.6 meV) from the 54 atom cell to the 432 cell and by 3.3 mRy 
(44.9 meV) from the 250 to the 432 atom cell. These changes in the doping energies as 



the 45 meV (3.3 mRy) 
6|, which suggests that 



the lattice size is increased are of the same order of magnitude as 
binding energy of the dopant as calculated by effective mass theory 
the effects of lattice relaxation and changes in the electronic structure will be important 
for the donor electron dynamics in phosphorus doped silicon where the isotropic doping 
fraction is ~ 0.2 % or higher. In systems with a lower doping fraction, these effects may 
still play a role, but we cannot make a definitive statement on this issue since the changes 



in energy exhibit a strong dependence on the system size up to the largest system studied 



here. In another recent study 32], the energy gained by relaxing a cell with a monolayer of 
phosphorus dopants was found to be of a similar magnitude. Whether the lattice relaxation 
is important to EDMR readout schemes depends on how the relaxation affects the scattering 
dynamics of electrons, which is related to the doping potential. In the next section, we make 
comparisons of the doping potentials for a donor phosphorus atom in silicon, both with and 
without the effects of geometrical relaxation of the lattice. 

A. Doping potential 

We define the doping potential as the doped cell potential minus the bulk silicon potential 
(Eqs. ([1]) and ([2])). The doping potential shows how the electronic environment surround- 
ing the dopant differs from that of bulk silicon. These potentials also provide input for 



calculations of the cross sections of electron scattering at the dopants 50|, l5l|, which are 



largely determined by integrals of the doping density [50|. By calculating the scattering of 
conduction electrons confined in a two-dimensional layer located at a given distance from 
the (001) plane, a connection can be made with electrically detected magnetic resonance 



(EDMR) schemes [33|, |3J, |52|-|54|] used to measure the dopant spin state. Doping potentials 
calculated using atomistic DFT can also be used to parameterize new tight binding models, 
or effective single-electron models which more accurately reproduce the effects of the dopant 
electronic structure than standard effective mass models. 

The doping potentials are shown for the 54 and 432 atom cells in Figs. [1] and |2l In the 
(001) plane, the doping potentials for the 54 atom cell (left panels) can be seen to overlap. 
Thus, while there are areas in which the doping potential goes to zero, indicating a return to 
bulk silicon behavior, the dopants are largely connected by regions of potential greater than 
100 meV. In the 432 atom cell at the bulk geometry, the doping potentials nearly go to zero 
between the dopants. However, there are small areas of non-zero potential still connecting 
dopants, suggesting that the single dopant limit has still not been reached. These potential 
connections between dopants are slightly exaggerated at the relaxed geometry. In both cells, 
however, the potential region directly around the dopants is similar, with a region exceeding 8 
eV directly around the dopant and three 1 eV lobes about 120 degrees apart from each other. 
In addition to the lobes near the nucleus, there are additional lobes at an energy scale of -700 
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FIG. 1: (Color online) The potential difference (eV) in the (001) plane between doped and undoped 



silicon is shown for cell sizes of 54 (1(a) and 1(c) ) and 432 atoms (1(b) and 1(d) ). An area of 93 x 



55 Ais shown for in all figures, in which the 54 and 432 atom cells repeat about 44 and 11 times, 
respectively. For the 54 and 432 atom cells, parts of 51 and 14 dopants, respectively, can be seen 



in the figures. Doping potentials are shown for cells at the bulk geometry (1(a) and 1(b)) and 



optimized geometries (1(c) and 1(d)). The contour lines are drawn every 1 eV except between -2 



and 2 eV, where they are drawn every 100 meV. The color axis in these figures is between -1 and 
1 eV in order to highlight the effects at this energy scale, while the range of doping energies is 
between -9 and 12 eV, with the larger values occuring near the dopant nucleus (see Fig. |3|). 

meV which extend away from the nuclei in space. The trigonal symmetry results from taking 
a two-dimensional cut through the three dimensional structure onto the (001) plane. When 
viewed in three dimensions, these lobes can be seen to arise from the tetrahedral nature 
of the bonding in the silicon cell. Three-dimensional potential isosurfaces at -700, -600, 
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FIG. 2: (Color online) The potential difference (eV) in the (001) plane between doped silicon in 
its optimized geometry and bulk silicon for the 432 atom cell in the bulk geometry. We refer to 



this as the "relaxing and doping potential." In Fig. 2(a), 14 phosphorus dopants are visible, while 



Fig. |2(b)] shows a close-up of the region around one dopant. The contours and color axis in Fig. 2(a 



are the same as in Fig. [TJ The close up in Fig. 2(b) uses the same contours, but the color axis is 
between -2 and 2 eV. 

and -150 meV are shown in Fig. [31 The non-spherical nature of the doping potential will be 
important for calculation of electron-dopant scattering cross sections. Effective mass models 
of doped silicon assume a spherical Coulomb potential, while some tight binding models use 



a spherical, Coulomb-like doping potential 28 



29 



55 



56| and allow the surrounding atoms' 



electronic structure to adjust according to this potential. The anisotropic doping potentials 
and cell geometries calculated here could be used as alternative parameterizations for tight 
binding models, and they can also be used to calibrate the resulting potentials and densities 
calculated by tight binding models. 

Fig. m shows a closer view of the doping potential in the region of the dopant for the 432 
atom cell with optimized geometry. In the effective mass picture {g], the dopant wavefunction 
has s-orbital character. In Fig. HI the potential in the region of the dopant can be seen to 
oscillate as a function of the orientation. The oscillations are due to interactions with 
electrons in the shells below the valence shell. If these calculations were performed without 
using pseudopotentials, which reduce oscillations from core electrons and replace them with 
a smooth potential, the potential would most likely oscillate to an even greater degree. These 
oscillations, as well as those visible in the optimized geometries of Fig. [1] due to the silicon 
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(a) 



(b) 





(c) 



FIG. 3: (Color online) Doping potential isosurfaces are shown for -700 ( |3(a)D , -600 ( |3(b)D , and 



150 (3(c) ) meV. The isosurfaces are shown in slightly transparent blue, and the phosphorus donors 
are shown in red. The silicon atoms are not pictured in Figs. 3(a) and |3(b)| and they are shown 



for reference in Fig. 3(c) 



lattice distortions, represent qualitative differences between DFT and effective mass theory. 
In Ref. 16|, Carter et. al. use mixed pseudopotentials in order to estimate the potential 



as a function of the distance from a layer of dopants. In Fig. [5l we plot the doping potential 



9 



SI^,P lo^fnizHf geometiyi Doqing Polsryiial (sV) 



Si^j^P {oplimiied geometry) Doping Po'ljerflial (eV] 




ho 



30 32 34 3e 3B 40 
Diseance along -cT O 0> (Angslronis) 

(a) 



33i 34 34.5 35 35.5 38 
Di&la^co along -cl C> {Angslf oms) 

(b) 



FIG. 4: (Color online) The potential difference (eV) between doped and undoped silicon is shown 
for the 432 atom cell at the optimized geometry of the doped cell for the region close to a dopant. 
The doping potential can be seen to oscillate in this region as opposed to showing s-orbital character 
as predicted by effective mass theory. Additionally, there are a number of features which are not 
spherically symmetric. The contour lines are the same as in Fig. [H while the color axis has been 
expanded to the range of -9 to 12 eV. 



(doped minus bulk) as a function of distance from the (001) layer of dopants. In Fig. 5(a 



full two-dimensional cut through the potential is given in the plane perpindicular to the (001) 



plane, and in Figs. 5(b) and 5(c) the potential is shown in a slice of this plane which connects 



two dopants. In contrast to Ref. 16|, where no structure is evident, marked structure is 
seen in Fig. |5l In Figs. 4 and 7 of Ref. 16[, the mixed pseudopotential doping potentials are 



much smoother than in Fig. 5(b) especially in the region around the dopant. In Fig. 5(b) 



there is a significant amount of structure in the potential near the dopant itself. Minor 



effects of the silicon atoms in the next layer of the crystal are also evident in Fig. 5(c) when 



the effects of geometric relaxation are included. It is important to note that these results 
are for doping densities near the single dopant limit: a study of the effect of a (5-layer of 
dopants would require very large cells which would likely have thousands of atoms. 
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FIG. 5: (Color online) The doping potentials iVsiiz-LP — ^5i) as a function of distance from the 



doping layer. Fig. 5(a) gives a contour plot with the x-direction representing an intralayer direction 
and the y-direction representing an interlayer direction. The contours and color axis are the same 
as in Fig. [TJ Fig. |5(b)] shows the potential as a function of distance from the doping layer along the 
line between dopants. The inset shows an enlargement for the energy range -200 to 200 meV, which 



is most relevant to scattering by a two-dimensional electron gas. Fig. 5(c) shows the potential with 
geometric relaxation included. 
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B. Constrained spin calculations 



In the effective mass model 



m 



19 



24l |. the dopants are treated as effective hydrogenic 



one-electron systems, with spin due to the additional dopant electron. Density functional 
theory gives a more detailed treatment of the many-body problem, allowing dopant electrons 
to couple to core electrons on the same dopant, to electrons in the silicon atoms, and most 
importantly, also allowing dopant electrons to couple with each other at low densities. In 
order to compare the DFT results with the frequently used single-electron picture 

m 

, we calculated the spin densities (spin up density minus spin down density) of 
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each doped cell. 

Fig. in] shows spin densities for DFT calculations in which the total spin of the unit cells 
was constrained to be = | a.u. (one bohr magneton). In all of the previous sections, 
the calculations were unconstrained in the spin degree of freedom and converged to 5 = 0, 
with zero spin density everywhere. In addition to giving an indication of spin ordering 
due to the coupling between dopant atoms in the system with respect to cell size, the spin 
density also gives a qualitative picture of where the "additional" electron provided by the 
phosphorus dopant is located when the total magnetization of the system is constrained. 
In the 54 and 128 atom unit cells, the spin density clearly shows interactions between the 
dopants in the form of areas of high density between dopants. The 250 and 432 atom cells 
show comparatively less localization of spin density between the dopant atoms, but there 
are still areas of enhanced spin density connecting the dopants along lines in the (001) 
plane. These areas are fading for the 432 atom cell but not completely absent. Note also 
that the maximum spin density is decreasing from the 54 to the 250 atom cell, but remains 
approximately constant between the 250 and 432 atom cells. These results give an effective 
one-electron picture of how modulations in the spin density are affected by the distance 
between dopant atoms and provide insight on the behavior of the effective one-electron 
wavefunction in this system. 

An exchange coupling between donors can be estimated using density functional the- 
ory according to Eq. 43-[49|. The quantity J provides an estimation of the spin cou- 



plings between donor qubits which may be used to apply two-qubit gates. The exchange 
couplings calculated for the different size supercells are given in Table [TTl As the dopant 
density decreases, we see that the magnitude of the exchange coupling also decreases. In 
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(e) 

FIG. 6: (Color online) DFT calculations of spin densities {n-^ — n^) of the 54, 128, 250, and 432- 
atom unit cells in a 93 x 55 Abox, with the spin constrained to S" = ^ in each unit cell. A close 



up of the dopant area of the spin density for the 432-atom cell is also shown (6(e) ). 



Refs. 



17| and (2l|, the exchange couphng of a two-donor system is calculated for dopants 



at much greater distances than in this work, using a Heitler-London approximation with 
variable alignment of the longitudinal and transverse Bohr radius of the dopant relative 
to the inter-dopant direction, respectively. These works, together with Refs. 23| and js^], 
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TABLE II: The exchange parameter as a function of ceU size. 



Unit Cell Distance between dopants (A) Exchange coupling (meV) 
SisaP 11.62 -125.1 
Sii27P 15.49 -83.6 
Si249P 19.36 -64.4 
Si43iP 23.23 -49.8 

have concluded that oscillations in the exchange coupling would make it difficult to control 
a quantum information system which attempts to exploit this coupling. Although donor 
spacings studied here are small compared to the previous studies of a pair of phosphorus 

n " 



donor atoms in bulk Si in Refs. [17| and [18|, we may nonetheless make a comparison of 



our results with these in Fig. [71 For the three largest systems studied here, the behavior 
of the exchange coupling with respect to donor spacing (r) is fit well by a single exponen- 
tial decay, J(r) = — 323.0 exp(—r/12.0). The exchange coupling at the distances studied 
using the Heitler-London approximation is systematically larger than the corresponding 
value extrapolated from our fitted decay. The source of the decrease in magnitude of the 
DFT estimated exchange coupling is likely due to oscillations in the donor electron densities 
which are present in the DFT calculation, but which are not included in the models used 



in Refs. [17[ and [18[. We note that due to the relatively small spacing between donors 
and the isotropic distribution of donors in this work, the correlations between donor atoms 
are stronger than what is modeled in the studies of isolated pairs of atoms separated by 
a large distance. However, if in actual devices the donor atoms are not evenly distributed 
and well separated, the magnitude of the exchange coupling may also be decreased with 
respect to what is predicted by the Heitler-London model. Extending our studies to larger 
systems with anisotropic doping will allow for a more direct comparison with studies based 
of isolated donor pairs. 

C. Doping density 



In Ref. 15[, wavefunctions of phosphorus dopants were calculated using effective mass 



theory, in which the Bloch functions of silicon where taken directly from a density func- 
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FIG. 7: Magnitude of the zero-field exchange couphng calculated by DFT in this work (circles) 



la] (tri- 



and in the Heitler-London approximation with effective mass theory wavefunctions in Ref. 
angles., Bohr radius of 2.381 nm), and in Ref. |l8| (squares, Bohr radius of 1.368 nm). Data for the 
three largest systems studied here are fit well by a single exponential decay in the spacing between 
donors (solid line.) The magnitude of the exchange coupling predicted by the fitted decay function 
for isotropic doping is much lower than that calculated in Refs. or 18] using the Heitler-London 
approximation for an isolated pair of P atoms in bulk Si (inset.) 

tional theory calculations. Doping densities (density of doped cell minus undoped cell) were 



presented in Fig. 2 of Ref 



15 



In Fig. IHl we present the doping density, 

P doping P doped Pundopedy 



(4) 



calculated for the 432 atom cell using a full atomistic DFT treatment, where each density p 
is calculated from the Kohn-Sham orbitals (pi as 



P 



(5) 



The current DFT doping densities show a different distribution in the (001) plane than 
those presented in Ref. 15|. In particular, the present densities are more circular around 
the dopant in this plane. The oscillations which are not included in the effective mass 
calculations are also evident in the immediate vicinity of the dopant. These oscillations are 
about two orders of magnitude less than the doping density in the vicinity of the dopants. 
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FIG. 8: (Color online) The doping density {psi^-nP — Psi) in the (001) plane. The region of the 



density near the dopant is zoomed in panel 8(b) The contours are drawn every 0.01 units, except 



between -0.1 and 0.0 where they are drawn every 0.002 units. The color axis is set between -6 and 
7 X 10~^ units, while the density varies between -6 x 10~^ and 0.35 units with the larger values 
near the nucleus. 



Lobes similar to those seen in the doping potentials in Section UlIAl are also apparent in the 
doping densities. Finally, we note that the doping density shows evidence of the interaction 
between dopants, resulting in finite electron density between dopants. 



IV. CONCLUSION 



We have presented density functional theory calculations for silicon doped with a single 
phosphorus atom, for systems with up to 432 atoms in a cell. A detailed knowledge of the 
electronic structure of doped silicon is necessary for the implementation of spin-based qubits 
in silicon [l-S, 5|. We have calculated non-spherical electron densities and doping potentials 



which allow for a level of microscopic description beyond effective mass and ti ght 
theory. In comparison to previous calculations of dopant electronic structure {gI-IigI 



28 



30 



bmdm 



17 



19 



3l|, we have found an unprecedented level of structure in the doping potentials 



(Fig. [T]) and densities/wavef unctions (Fig. |8]). The exchange coupling between qubits was 



16 



found to be less than estimates based on the Heitler-London approximation, due to the 
oscillatory nature of doping potentials. 

These calculations have many potential uses. Such detailed microscopic calculations 
will allow more accurate and detailed device simulations than are currently possible. By 
understanding the effects of modulations in the doping density including effects of both 
the spin density as well as the doping potential, allows now calculations which probe the 
readout properties of the qubits, especially using techniques such as EDMR 



33 



34 



52 



54| . Additionally, alternative qubits such as excited-state dopants |58| or charged dopant 



qubits 



56 



59 



|-|6lj can now be explored with this accurate picture of the electronic structure. 
The doping potentials calculated here can provide the starting point for effective one-electron 
calculations of a dopant electron wavefunction, possibly deformed by some electrostatic 
gate potential. They also can guide design of multiple qubit devices by providing effective 
Hamiltonians or potentials for multiple-qubit geometries. Finally, the doping potentials 
provide input for scattering calculations, including calculations in which the current-carrying 
electrons are confined to a two-dimensional plane to model electrical readout schemes for 
silicon quantum computation. The spin densities we have calculated here can also be used 
to compare with a single-electron picture and to determine the density of the electron which 
is donated by the phosphorus. 

In the future, we plan to look at systems which more accurately represent the experimental 
devices. This will require looking at the effect of the silicon dioxide interface and defects 
at this interface, a considerably more computationally intensive task. They may also be 
extended to the calculation of parameters related to the hyperfine interaction. These DFT 
calculations may also be coupled with calculations of the spin-dependent scattering 

1 as quantum control calculations for the 



5l| of the two-dimensional electron gas as we 
implementation of quantum logic operations 



62 



66|. 
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